The impact of phage and phage resistance on microbial community dynamics

Where there are bacteria, there will be bacteriophages. These viruses are known to be important players in shaping the wider microbial community in which they are embedded, with potential implications for human health. On the other hand, bacteria possess a range of distinct immune mechanisms that provide protection against bacteriophages, including the mutation or complete loss of the phage receptor, and CRISPR-Cas adaptive immunity. While our previous work showed how a microbial community may impact phage resistance evolution, little is known about the inverse, namely how interactions between phages and these different phage resistance mechanisms affect the wider microbial community in which they are embedded. Here, we conducted a 10-day, fully factorial evolution experiment to examine how phage impact the structure and dynamics of an artificial four-species bacterial community that includes either Pseudomonas aeruginosa wild-type or an isogenic mutant unable to evolve phage resistance through CRISPR-Cas. Additionally, we used mathematical modelling to explore the ecological interactions underlying full community behaviour, as well as to identify general principles governing the impacts of phage on community dynamics. Our results show that the microbial community structure is drastically altered by the addition of phage, with Acinetobacter baumannii becoming the dominant species and P. aeruginosa being driven nearly extinct, whereas P. aeruginosa outcompetes the other species in the absence of phage. Moreover, we find that a P. aeruginosa strain with the ability to evolve CRISPR-based resistance generally does better when in the presence of A. baumannii, but that this benefit is largely lost over time as phage is driven extinct. Finally, we show that pairwise data alone is insufficient when modelling our microbial community, both with and without phage, highlighting the importance of higher order interactions in governing multispecies dynamics in complex communities. Combined, our data clearly illustrate how phage targeting a dominant species allows for the competitive release of the strongest competitor while also contributing to community diversity maintenance and potentially preventing the reinvasion of the target species, and underline the importance of mapping community composition before therapeutically applying phage.


Introduction
Microbiome research is a dynamic and growing field in microbiology, producing an incredible amount of sequence data from a wide range of clinical and environmental samples.Humans, for instance, are colonised by a large number of microorganisms and research continues to implicate microbial communities as potential drivers behind multiple important biological processes [1][2][3].These processes may play important roles in human health and disease, with some work focusing on correlations based on microbiome composition [4][5][6][7][8] while other look more closely for direct causality [9][10][11][12].Still, the challenge to move beyond descriptive and correlative approaches remains, and there is a need to develop bottom-up mechanistic and quantitative understanding of the forces acting upon and shaping microbial communities.To this end, synthetic polymicrobial communities are being designed, and are gaining traction in both pure and applied microbiome studies [13][14][15][16].Synthetic microbiomes allow for precise and malleable experimental testing of the basic rules that govern both microbial organisation and functioning on molecular and ecological scales [17][18][19][20], as well as allowing for exploration of causal roles connecting specific microbiome structures to potential outcomes of interest.
Bacteria and their viral predators, bacteriophages (phages), have long been of interest in microbiological research, in part due to being the most abundant biological entity on the planet [21,22].Phages are highly diverse in terms of their morphology, genetics, and life histories [21,23], with a clear distinction between obligatory killing lytic phages and temperate phages that can either cause a dormant infection (lysogenic cycle) or cell lysis to release new phage particles (lytic cycle).Phages are thought to play a key role in shaping both the taxonomic and functional composition of microbial communities as well as their stability, ecology, and evolution [23][24][25][26][27].For example, lytic replication will per definition cause a reduction in the density of the bacterial host strain or species, which in turn can have knock-on effects for the microbial community composition through the enabling of invasion and/or coexistence of competitor species.Despite the large potential impact of lytic phage, only a very limited number of experimental studies have explored the ecology and evolution of bacteria-phage interactions in a microbial community context [28,29], and it remains unclear if and how interactions between different species in more complex communities shape the effects of lytic phages on microbial eco-evolutionary dynamics.Consequently, we lack the stepping stones to understand how phages shape microbial community dynamics (reviewed in [23]), which are urgently needed to understand potentially causal relationships between natural phage communities and a variety of human diseases [30][31][32][33][34][35], and for optimising the therapeutic application of phages in the clinic.
A key consideration in this context is that bacteria can overcome phage infection through a range of different means [36,37], with varied underlying molecular mechanisms and which can act during different stages of phage infection [38][39][40][41].Through the modification, masking or complete loss of phage-binding surface receptors, for example, bacteria can prevent phage adsorption and injection [40,42].Systems such as CRISPR-Cas on the other hand work by inserting short DNA sequences from phage and other invasive mobile genetic elements into the host genome to provide future immunological memory [43].Unlike CRISPR-based resistance [15], phage resistance through receptor mutation can be associated with substantial phenotypic shifts and fitness trade-offs, through changes to virulence [44,45], biofilm formation [46], or antibiotic resistance [47].
While our previous work asked how interspecific competition shapes phage resistance evolution in Pseudomonas aeruginosa [15], we here sought to answer the inverse and complimentary question of how host-phage interactions shape the composition and stability of the wider microbial community.To this end, we combined exploratory and hypothesis-driven approaches, applying experimental evolution to examine how a phage impacts the dynamics of an artificial bacterial community.This community consisted of Pseudomonas aeruginosa, Staphylococcus aureus, Acinetobacter baumannii, and Burkholderia cenocepacia, all of which are opportunistic pathogens that can cause severe infection and may coinfect with one another [48][49][50][51].Firstly, we hypothesised that the addition of a P. aeruginosa-specific phage would promote species coexistence by limiting P. aeruginosa dominance through competitive release (expansion of phage-resistant competitors, following removal of phage susceptible competitor) in a way akin to what is commonly observed with antibiotics [14,[52][53][54][55]. Secondly, we hypothesised that blocking the ability of P. aeruginosa to evolve CRISPR-based immunity would reduce P. aeruginosa persistence due to community-dependent fitness costs of surfacemodification [15].Finally, we asked if we could quantitatively capture community dynamics using mathematical modelling.We found that the addition of a P. aeruginosa targeting phage resulted in the general maintenance of community diversity and coexistence, but also a shift in dominant species from P. aeruginosa to A. baumannii-with the former being unable to reinvade even after the phage was driven extinct.The impact of the type of phage resistance was limited or transient; however, while a P. aeruginosa wild type (WT) with the ability to evolve CRISPR-based phage resistance initially had a slight fitness advantage in the presence of A. baumannii over its CRISPR-negative isogenic mutant, this effect was lost over time as the phage was driven extinct.

Results
To measure the effect of phage on microbial community dynamics, we carried out a fully factorial 10-day in vitro evolution experiment using all possible combinations of 1, 2, 3, or 4 competitor species: S. aureus, A. baumannii, B. cenocepacia, and P. aeruginosa PA14 in the presence or absence of lytic phage DMS3vir.We previously applied the same model community to explore the effect of interspecific competition on phage resistance evolution in the P. aeruginosa WT over 3 days in the presence of phage [15].Here, we include both the WT P. aeruginosa PA14 strain, which can evolve CRISPR-based phage resistance, and an isogenic mutant lacking a functional CRISPR system to examine the impact of CRISPR-Cas versus surface modification on these dynamics.Following inoculation, we tracked the microbial community dynamics for all experimental treatments at regular intervals over a period of 10 days.All experiments were conducted in lysogeny broth (LB) at 37˚C (see Methods for details).

P. aeruginosa dominates in the absence of phage
Without P. aeruginosa present in the community, S. aureus was primarily the dominant species-with the ability to coexist with A. baumannii while outcompeting B. cenocepacia (Figs 1  and S1).This, however, was not reflected once P. aeruginosa was introduced to the community.In the absence of phage, P. aeruginosa quickly became the dominant species in the microbial community, regardless of starting composition and the P. aeruginosa genotype (PA14 WT versus CRISPR-KO) (Figs 2 and S3).Consistent with this, the densities of the competitor species rapidly declined during these coculture experiments (Fig 2).Yet, there was a clear difference in the rate at which competitor species declined in frequency, which was highest for S. While the microbial community dynamics were relatively similar for the WT and CRISPR-KO strains, some significant differences were observed.For example, the densities of the CRISPR-KO strain were slightly lower in the presence compared to the absence of S. aureus on its own (Fig

Phage affects microbial community dynamics
Whereas P. aeruginosa dominated in the absence of phage, we hypothesised this would change once a PA14 targeting phage (DMS3vir) was introduced, largely by a virulent phage reducing the susceptible host population, facilitating expansion of other species through competitive release [14,[52][53][54][55].As expected, phage DMS3vir initially reached high titres due to replication on sensitive P. aeruginosa hosts, followed by a rapid decline in phage densities due to the evolution of phage resistance, regardless of whether the host had a functional CRISPR-Cas system or not (Fig 3C  Interestingly, the PA14 WT generally reached greater relative abundance than the CRISPR-KO strain when in the presence of A. baumannii, consistently doing so early in the experiment when phage remained in the population (Figs 3A, 3B, and S2).This was in Total abundance (cfu/mL)  concordance with P. aeruginosa evolving higher levels of CRISPR-based immunity against phage DMS3vir in treatments including A. baumannii due to the increased fitness cost of surface modification (Fig 4 and [15]): At 3 days postinfection, there was a significant effect of all treatments on the proportion of CRISPR-based resistance that had evolved compared to the PA14 monoculture, but this effect was strongest for treatments that contained A. baumannii.
At time point 10 we only found an increased proportion of P. aeruginosa clones immune through CRISPR-Cas when the treatment included A. baumannii (GLM; A. baumannii; t = 2.637, p = 0.01; S. aureus and A. baumannii, t = 2.283, p = 0.025; A. baumannii and B. cenocepacia, t = 2.689, p = 0.0087; polyculture, t = 2.141, p = 0.035).Overall, however, it was evident that mutation of the Type IV pilus became the dominant resistance mechanism even if P. aeruginosa has a functional CRISPR system (Fig 4), which might in part be why P. aeruginosa did not recover in the microbial community post phage exposure due to the associated fitness costs [15].

The type of evolved phage resistance does not have a knock-on effect on microbial community dynamics
Phage DMS3vir targets P. aeruginosa's Type IV pilus (T4P), an important virulence factor [56].We have previously shown that the evolution of phage resistance by mutation of the pilus is associated with large fitness trade-offs in the same microbial community as used in this study, whereas evolution of CRISPR-based immunity is not associated with any detectable trade-offs [15].We therefore predicted that the ability to evolve phage resistance through CRISPR-Cas would also have knock-on effects for the microbial community dynamics.However, measurement of the abundance of the competitors revealed that these were overall largely unaffected by the presence of a functional CRISPR-Cas immune system in P. aeruginosa with the exception of S. aureus:

A P. aeruginosa targeting phage results in the competitive release of A. baumannii and general diversity maintenance
We hypothesised that the effect of phage on microbial community structure could largely be explained by the competitive release (increase in absolute abundance, following removal of competitor) of A. baumannii, which then takes over to become the dominant species [57].To assess this, we examined the fold change difference for the final abundance of all 3 community members in the presence versus absence of phage (Fig 5).Crucially, this revealed a strong increase in A. baumannii density in the presence of a phage, supporting the idea that it becomes the dominant and determinant community member when P. aeruginosa is inhibited by phage (Fig 3).By contrast, when phage was added, S. aureus only experienced a clear fold change increase if it was cocultured with the CRISPR-KO strain and an additional competitor species.B. cenocepacia meanwhile seemed to be the species with the least benefit of phage, but still with a small fold change increase for some treatments (Fig 5).
The substantial fold increase in A. baumannii given the presence of phage (Fig 5) reflects a sustained divergence in the trajectory of A. baumannii in the phage treatments, despite the attenuation of phage titre by day 7 (Fig 3C).We hypothesised that the lack of P. aeruginosa rebound after phage clearance was due to a frequency-dependent shift in competitive dominance.To test this hypothesis, we competed ancestral A. baumannii, S. aureus, and B. cenocepacia against increasingly rare P. aeruginosa challenge, and found no barrier to P. aeruginosa invasion in pairwise experiments, down to a frequency of 1 in 10,000 cells (Fig 6).This result suggests that the failure of P. aeruginosa to return to dominance following phage clearance is due to more complex community-mediated interactions.
Additionally, we tested if it was A. baumannii that could have gained an advantage through natural selection when competing against P. aeruginosa over time, with phage allowing A. baumannii to better adapt to the environment and explain the inability of PA14 to reinvade.Yet, our data demonstrated that there was no difference in competitive performance of evolved A. baumannii relative to its ancestral strain, and both the ancestor and the evolved clonal populations of A. baumannii were outcompeted by the PA14 wild type ( Fig 7).
To quantitatively assess changes in community diversity, we calculated Shannon diversity indexes for all experimental treatments.We hypothesised that the addition of phage not only results in competitive release of one other bacterium (Fig 5), but also facilitates general maintenance of microbial diversity.Plotting these diversity scores over time shows that without phage, there is a rapid loss of diversity over time, whereas community complexity persists in the presence of phage ( ).This was true for treatments for both P. aeruginosa genotypes, but the trend became most pronounced for the CRISPR-KO strain when applying direct comparisons using Tukey contrasts, in which case we found phage to significantly increase diversity over time in nearly all treatments (

Four-species community dynamics are predictable from 2 and 3-species community data, in the absence of phage
Mathematical modelling provides a platform to quantify ecological interactions that determine community-level behaviours, as well as identify rules governing qualitative system behaviour [60][61][62][63].However, a major challenge in synthetic community research is developing robust modelling frameworks that are capable of predicting community dynamics [64].In a final set of analyses, we sought to parameterise and assess the predictive performance of generalised Lotka-Volterra (gLV) competition equations, trained on just 2-species data or a combination of 2-and 3-species data.Our results showed that fitting gLV models with pairwise only datasets led to predictive failures when applied to 3-or 4-species datasets (S6 Fig), consistent with the presence of higher order interactions effects (when the effect of species A on species B is dependent on the presence of species C [65]).In contrast, fitting gLV models to 2-and 3-species data and using the resulting interaction terms to predict 4-species dynamics reasonably fit the data in the absence of phage (Fig 9; fitted model coefficients are in S8 Fig).
In the presence of phage (Fig 10), we again utilised the gLV framework where the impact of phage is implicit (quantified by how interaction coefficients change as compared to the nophage case).The gLV model framework could adequately describe 2-and 3-species data, but the interaction coefficients did not generalise quantitatively to 4-species data-likely reflecting the structural limitation of a gLV competition model that does not explicitly capture phage predation dynamics.However, the model parameterised with 1-, 2-, and 3-species data did capture a qualitative shift in ecological outcomes from sole P. aeruginosa survival to competitive release of A. baumannii and S. aureus when P. aeruginosa is targeted by phage (S7 Fig).

Discussion
The advent of deep sequencing has dramatically increased our knowledge of the composition and functioning of microbiomes both in and around us.The role of microbial communities in human health has consequentially received increasing attention, with research focusing on how changes in microbiome composition over time may affect human health and define patient outcomes (reviewed in [66]).In addition, an increasing number of correlational studies find associations between virome composition and the health status of their host [23,[67][68][69][70], likely mediated by changes in the microbiome that could be either cause or effect.A deeper understanding of the impact of phages on microbiomes is likely to help to infer causal relationships between viromes and human health, and to design optimal therapeutic phage interventions (phage therapy).
Here, we expanded on our previous work on how interspecific competition can shape the evolution of phage resistance in a focal species (P.aeruginosa) [15], to study how the interaction between phage and bacterial immune mechanisms affects the broader microbial community dynamics.We found that whereas P. aeruginosa dominated in the absence of a phage, the presence of phage resulted in microbial diversity maintenance and A. baumannii becoming the dominant species (Figs 3 and 5).Interestingly, the competitive release of A. baumannii occurred in all treatments and was virtually independent of whether P. aeruginosa had a functional CRISPR-Cas immune system or not.This showed that the amplification of the fitness cost of P. aeruginosa receptor mutation in the presence of competitor species [15] has limited impact on the overall community dynamics.Overall, our experimental data align with the notion of phages having the potential to increase diversity and microbiome stability [27,71,72], and support the idea that phages can be useful in the designing of synthetic microbial communities [73].Surprisingly, our data do not support the hypothesis that bacterial adaptive immune systems play an important role in phage-mediated microbial community structuring under the experimental conditions tested here.
Our mathematical analyses focused on the ability of gLV models to predict community dynamics.While our analyses showed reasonable predictive success when incorporating Model for no phage data.Model fit predictions for 2-, 3-, and full 4-species community dynamics (solid lines) compared to experimental data (dashed lines) (PA = P. aeruginosa, SA = S. aureus, AB = A. baumannii, BC = B. cenocepacia).Models of 2-and 3-species dynamics were parameterised via optimization with least-squares to fit to a system of ODEs (defined as a gLV competition model with n species, where n = 1, 2, 3, 4).Only single species maximal growth rates (r i for species i = 1,. ..,n) were fixed from fitting mono-culture data, all interaction coefficients (β i,j describing the inhibitory effect of species j on species i for all i,j = 1,2 in the 2-species case, for all i,j = 1,2,3 in the 3-species case) were open for fitting.We construct the full 4-species community interaction matrices (one for PA14-shown in S8 Fig -and one for CRISPR-KO) by averaging corresponding β i,j interaction terms from the fit 2-and 3-species models (see S1 Text), and use this matrix to simulate dynamics in the respective polyculture cases.See Methods and S1 Text for detailed description of mathematical modelling.The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.https://doi.org/10.1371/journal.pbio.3002346.g009Here, models were parameterised via optimization with least-squares to fit a system of ODEs (defined as a gLV competition model with n species, where n = 1, 2, 3, 4), where we do not explicitly track the phage population dynamics.Instead, we assume that the phage acts as some external perturbation that leads to changes in the interactions between community members (β i,j values differ from values in Fig 9).Only single species maximal growth rates were fixed from fitting mono-culture data (r i for species i = 1,. ..,n), all interaction coefficients were open for parameterising 2-and 3-species models (β i,j describing the inhibitory effect of species j on species i for all i,j = 1,2 for the 2-species case, for all i,j = 1,2,3 for the 3-species case).We then construct the full 4-species community interaction matrices (one for PA14 and one for CRISPR-KO) by averaging corresponding β i,j interaction terms from the fit 3-species models (treatments: PA+AB+SA, PA +BC+SA, PA+AB+BC with phage), and use this matrix to simulate dynamics in the respective polyculture cases.See Methods and S1 Text for detailed description of mathematical modelling.The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.https://doi.org/10.1371/journal.pbio.3002346.g0103-species data, we note that our analyses pose 2 distinct questions: (1) How can we provide more accurate predictions?(2) What general lessons can we draw from our model analyses?
In agreement with a growing number of gLV-based analyses, we found that a simple "bottom up" model fitting approach (fitting single species growth, then all pairwise interactions, then predicting larger system behaviour [14]) performed poorly, indicating the presence of significant higher order interactions [65,74,75].Consistent with this conclusion, we found that allowing pairwise interactions to vary (contingent on the presence of a third species) produced both qualitative and quantitative improvements in predicting community dynamics (Figs 9 and 10).In the presence of phage, our model successfully predicted the qualitative result of A. baumannii competitive release, but failed to quantitatively replicate observed community dynamics (Fig 10).This quantitative failure suggests that our underlying gLV model structure (the dominant framework in microbiome modelling studies [61,71,76]) excludes critical components, such as higher order and/or heterogeneous (in time or space) interactions as well as the explicit predatory effect of phage on P. aeruginosa (also likely time and spatially dependent).Additionally, it emphasises an ongoing need in microbiome modelling to evaluate functional forms that can efficiently-with respect to parameter number-and accurately capture the complexities of community dynamics.
Our parameterised models are tuned to the data generated by our specific 4-species community, which raises the question of "can we learn more general lessons from our model?"If we simplify our analysis to a 2-species context (focal pathogen, subject to phage, plus a second, non-focal species), we can translate recent analyses on the impact of (antibiotic) perturbations in a 2 species context [55].This approach delivers a couple of general messages.First, we can provide a general mathematical definition of "competitive release" mediated by phage predation (see S1 Text) highlighting the importance of both demographic and species interaction parameters.Second, we can underline that phage control of a focal pathogen presents secondary ecological problems, if the pathogen is competing with other pathogens that are not targeted by the phage.In this scenario, phage therapy (or other "narrow spectrum" treatment) can lead to competitive release of previously rare pathogens, as seen in our experimental data showing the replacement of P. aeruginosa by A. baumannii, following phage treatment.These results imply that "narrow spectrum" anti-microbials, such as phages, may not always be the best option when multiple pathogen species are competing within a single polymicrobial infection.One counter-intuitive suggestion, grounded in the idea of "beneficial resistance" [55], is to co-administer probiotic competitors that are resistant to the treatment (i.e., phage or antibiotic resistant) and can therefore continue to exert ecological suppression on the focal pathogen during the course of treatment, while presenting minimal direct risk of disease.Alternatively, one could apply phage cocktails that target not just the dominant pathogen, but also other coexisting bacterial pathogens, to preemptively prevent their invasion.

Bacteria and phages
The bacteria P. aeruginosa UCBPP-PA14 strain marked with streptomycin resistance, the PA14 csy3::LacZ strain (CRISPR-KO), and phages DMS3vir and DMS3vir+acrF1 were used throughout this study and have all been previously described [77,78].The microbial community consisted of S. aureus strain 13 S44 S9 and A. baumannii clinical isolate FZ21 which were isolated at Queen Astrid Military Hospital, Brussels, Belgium, while B. cenocepacia J2315 was originally isolated from a person with cystic fibrosis in the UK in 1989 and was provided by Queen Astrid Military Hospital, Brussels, Belgium.
(Cetrimide agar (Invitrogen) for PA14 selection and LB agar supplemented with 50 μg/ml of gentamicin to select for A. baumannii) and 6 randomly selected colonies from each replicate of indicated treatment/time point were pooled and inoculated overnight in 6 ml of LB medium at 37˚C with agitation (n = 6 per treatment, unless indicated otherwise).In parallel, 6 colonies from the ancestral strains were pooled and subject to the same overnight growth conditions.After 24 h of growth, competition assay and sample treatments were performed as described above.
To determine the competitive performance of the focal species relative to competitor strain, we used the selection rate (r), defined as the difference in Malthusian parameters as follows: r = (ln[density strain A at day t/density strain A at t−1]-ln[density strain B at day t/density strain B at t−1])/day) [58,59].The data used for these calculations were the bacterial quantities (cfu/ml) as estimated by qPCR as explained above, with 2 technical replicates per assay.

Mathematical modelling
Models were parameterised via optimization with least-squares regression to fit the gLV competition model, , where N i (t) is the density of the ith species, r i is the respective single species maximal growth rate, β ij describes the per capita effect of species j on species i, and n is the total number of species.We take a "bottom up" approach [81] to determine the interaction coefficients β i,j .In all cases, we determine single species maximal growth rates r i from mono-culture time series data and fix them for 2-, 3-, and 4-species model parameterisation.Initially, we fit pairwise interaction coefficients for all possible 2-species cocultures and from here, construct an interaction matrix to predict the dynamics for the 3-and 4-species communities (S6 Fig).This is done for both PA14 and CRISPR-KO strains, with and without phage, where phage effects are implicitly represented by changes in interaction parameters between the models with and without phage.To improve results, we additionally fit pairwise interaction parameters β i,j using 3-species experimental data where all interaction parameters are open (only growth rates fixed).Using either the resulting interaction terms or averaging these coefficients with the 2-species coefficients (in PA14 no phage case, S1 Text), we are again able to construct an interaction matrix to predict 4-species community dynamics (Figs 9 and 10).
See S1 Text for further description of above model parameterisation methods, simulation methods (S7 Fig), and mathematical analysis of phage-dependent competitive release.All modelling and analysis was done using Matlab 2021b and the code is publicly available at: https://github.com/GaTechBrownLab/phage-community-dynamics.git.

Statistical analyses
Analysis of the effects of the various species compositions on P. aeruginosa densities in the absence (Fig 2 ) or presence (Fig 3) of phage were done using a generalised linear model (GLM) approach, with log10 cfu/ml set as the response variable.The explanatory variables used in the analyses were type of PA14 clone (PA14 WT or CRISPR-KO), treatment, time point, replica, and experimental repeat to account for potential pseudo-replication.
To explore the impact of interspecific competition on the evolution of phage resistance at time points 3 and 10 (Fig 4), we used a quasibinomial GLM where the proportion of evolved CRISPR-based phage resistance was the response variable, and treatment, replica, and experimental repeat were the explanatory variables.
The analyses of fold-changes to assess competitive release by comparing absolute density differences of the individual community members in the absence v presence of phage (Fig 5 ; S. aureus, A. baumannii, and B. cenocepacia) was done through Wilcox signed rank exact tests.A nonparametric test was chosen after performing a Shapiro-Wilk test for normality.
Next, the diversity maintaining effects were examined through assessing the effect of phage DMS3vir on Shannon Diversity index scores over time (Fig 8).This was done through a linear model where the Shannon Diversity index score (H) was the response variable, and treatment, time point, the presence of phage, PA14 clone (PA14WT and CRISPR-KO), experimental repeat, and replica were the explanatory variables.Shannon Diversity (H), was calculated as H = -Sp i * ln(p i ), where S is the sum and p i is the proportion of the entire community made up of species i.
For the competition assay (Figs 6 and 7), Graphpad Prism9 software (San Diego, California, USA) was used for statistical analysis.We used one-way ANOVA with Tukey post hoc testing for multiple comparisons, in which p < 0.05 was considered statistically significant.
Throughout the paper, pairwise comparisons were done using the Emmeans package [82], and model fits were assessed using Chi-squared tests and by comparing Akaike information criterion (AIC) values, as well as plotting residuals and probability distributions using histograms and quantile-quantile plots (Q-Q plots), respectively.All statistical analyses were done using R version 4.3.0.[83], its built-in methods, and the Tidyverse package version 2.0.0 [84].All data is available at: https://doi.org/10.6084/m9.figshare.24187284.cenocepacia).Models were parameterised via optimization with least-squares to fit a system of ODEs (defined as a generalised Lotka-Volterra competition model with n species, where n = 1, 2, 3, 4).We parameterise the models via fitting of 1-(for growth rates r i ) and 2-(for all possible pairwise interaction coefficients β i,j 8i,j = 1,2) species dynamics and use the resulting coefficients to predict the 3-and 4-species community dynamics.For fitting coculture data, growth rates r i were fixed from mono-culture data and interaction parameters β i,j were all open.See Methods and S1 Text for a detailed description of mathematical modelling.The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.
For visualisation purposes, the data from Fig 3 is also presented as an ordination plot (S5 Fig).

Fig 2 .
Fig 2. P. aeruginosa becomes the dominant species in the absence of phage.Showing the community composition and bacterial densities in cfu/ml over time for the microbial communities in the absence of phage for the communities with either the (A) PA14 WT or (B) CRISPR-KO mutant.The community composition was estimated by qPCR at days 0, 1, 3, 7, and 10 of the experiment.The coloured bars represent the relative abundance of each species (left y axis), while the white line represents total abundance in cfu/ml (right y axis).Each panel represents average composition across 6 replicates for each treatment over time.PA14 = P. aeruginosa, SA = S. aureus, AB = A. baumannii, BC = B. cenocepacia, MC = microbial community.For individual replicates of species abundance, see S3 Fig.The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.https://doi.org/10.1371/journal.pbio.3002346.g002

Fig 3 .
Fig 3. Phage allows for the maintenance of all microbial community members, with A. baumannii becoming the new dominant species.Showing the community composition and bacterial densities in cfu/ml over time for the microbial communities in the absence of phage for the communities with either the (A) PA14 WT or (B) CRISPR-KO mutant.For A and B, the community composition was estimated by qPCR at days 0, 1, 3, 7, and 10 of the experiment.The coloured bars represent the relative abundance of each species (left y axis), while the white line represents total abundance in cfu/ml (right y axis).Each panel represents average composition across 6 replicates for each treatment over time.PA14 = P. aeruginosa, SA = S. aureus, AB = A. baumannii, BC = B. cenocepacia, MC = microbial community.For individual replicates of species abundance, see S4 Fig. (C) Phage titres for phage DMS3vir over time across all experimental treatments (PA = P. aeruginosa, SA = S. aureus, AB = A. baumannii, BC = B. cenocepacia, MC = microbial community), infecting either the PA14 WT or the CRISPR-KO strain as indicated by line type.Each data point represents a replicate, with lines following the mean and the error bars denoting 95% CI.Asterisks indicate a significant overall difference in phage density between the PA14 WT (n = 12 per time point) or CRISPR-KO clone (n = 6 per time point) (effect of P. aeruginosa clone; linear model: * p < 0.05).The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.https://doi.org/10.1371/journal.pbio.3002346.g003

Fig 5 .
Fig 5. Fold change between no phage and phage treatments at the end of the experiment.The fold change difference of the individual community species not targeted by phage when comparing absolute densities in the presence of phage to the absence at the final experimental time point.Asterisks indicate higher final absolute density in the presence versus absence of phage (Wilcoxon signed rank exact test: * p < 0.05, ** p < 0.01, *** p < 0.001).The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.

Fig 6 .Fig 7 .Fig 8 .
Fig 6.P. aeruginosa can invade from low initial frequency against all community members.Showing P. aeruginosa density in cfu/ml from competition experiments between PA14 wild type with variable starting densities against either S. aureus (SA), A. baumannii (AB), or B. cenocepacia (BC).The species densities were estimated by qPCR at time point 0 and 24 h post coculture.Box plots show the median, 25th and 75th percentile, and the interquartile range.Raw values from each biological replicate are shown as points (n = 6 per pairwise competition).The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.https://doi.org/10.1371/journal.pbio.3002346.g006

Fig 9 .
Fig 9.Model for no phage data.Model fit predictions for 2-, 3-, and full 4-species community dynamics (solid lines) compared to experimental data (dashed lines) (PA = P. aeruginosa, SA = S. aureus, AB = A. baumannii, BC = B. cenocepacia).Models of 2-and 3-species dynamics were parameterised via optimization with least-squares to fit to a system of ODEs (defined as a gLV competition model with n species, where n = 1, 2, 3, 4).Only single species maximal growth rates (r i for species i = 1,. ..,n) were fixed from fitting mono-culture data, all interaction coefficients (β i,j describing the inhibitory effect of species j on species i for all i,j = 1,2 in the 2-species case, for all i,j = 1,2,3 in the 3-species case) were open for fitting.We construct the full 4-species community interaction matrices (one for PA14-shown in S8 Fig -and one for CRISPR-KO) by averaging corresponding β i,j interaction terms from the fit 2-and 3-species models (see S1 Text), and use this matrix to simulate dynamics in the respective polyculture cases.See Methods and S1 Text for detailed description of mathematical modelling.The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.

Fig 10 .
Fig 10.Model for phage data.Model fit predictions for 2-, 3-, and full 4-species community dynamics (solid lines) in the presence of phage compared to experimental data (dashed lines) (PA = P. aeruginosa, SA = S. aureus, AB = A. baumannii, BC = B. cenocepacia).Here, models were parameterised via optimization with least-squares to fit a system of ODEs (defined as a gLV competition model with n species, where n = 1, 2, 3, 4), where we do not explicitly track the phage population dynamics.Instead, we assume that the phage acts as some external perturbation that leads to changes in the interactions between community members (β i,j values differ from values inFig 9).Only single species maximal growth rates were fixed from fitting mono-culture data (r i for species i = 1,. ..,n), all interaction coefficients were open for parameterising 2-and 3-species models (β i,j describing the inhibitory effect of species j on species i for all i,j = 1,2 for the 2-species case, for all i,j = 1,2,3 for the 3-species case).We then construct the full 4-species community interaction matrices (one for PA14 and one for CRISPR-KO) by averaging corresponding β i,j interaction terms from the fit 3-species models (treatments: PA+AB+SA, PA +BC+SA, PA+AB+BC with phage), and use this matrix to simulate dynamics in the respective polyculture cases.See Methods and S1 Text for detailed description of mathematical modelling.The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.

S1
Text.Description of the mathematical modelling.(PDF) S1 Fig. Line plot of bacterial densities in the absence of P. aeruginosa and its phage.Showing the bacterial densities in cfu/ml over time for SA (S. aureus), AB (A. baumannii), and BC (B.cenocepacia) in various coculture combinations in the absence of P. aeruginosa and its phage.Dashed horizontal line at 10 2 cfu/ml marks the threshold of reliable detection where the qPCR results indicate the bacteria has gone or is close to extinction from a population.Data are mean ± 95% CI.The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.(EPS) S2 Fig. Ordination plot in the absence of phage.PCA ordination of relative bacterial abundance in the absence of phage DMS3vir, with grid layouts separated into days post phage infection.Outer circle colour indicates which PA14 clone is present in the population, while inner circle indicates community composition (SA = S. aureus, AB = A. baumannii, BC = B. cenocepacia).The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.(EPS) S3 Fig. Line plots of bacterial densities in the absence of phage.Showing the bacterial densities in cfu/ml over time for the PA14 WT and CRISPR-KO P. aeruginosa strains, and b the other microbial community species (SA = S. aureus, AB = A. baumannii, BC = B. cenocepacia, MC = microbial community) in the absence of phage DMS3vir.Dashed horizontal line at 10 2 cfu/ml marks the threshold of reliable detection where the qPCR results indicate the bacteria has gone or is close to extinction from a population.Data are mean ± 95% CI.The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.(EPS) S4 Fig. Line plots of bacterial densities in the presence of phage.Showing the bacterial densities in cfu/ml over time for the PA14 WT and CRISPR-KO P. aeruginosa strains, and b the other microbial community species (SA = S. aureus, AB = A. baumannii, BC = B. cenocepacia, MC = Microbial community) in the presence of phage DMS3vir.Dashed horizontal line at 10 2 cfu/ml marks the threshold of reliable detection where the qPCR results indicate the bacteria has gone or is close to extinction from a population.Data are mean ± 95% CI.The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.(EPS) S5 Fig. Ordination plots in the presence of phage.PCA ordination of relative bacterial abundance in the presence of phage DMS3vir, with grid layouts separated into days post phage infection.Outer circle colour indicates which PA14 clone is present in the population, while inner circle indicates community composition (SA = S. aureus, AB = A. baumannii, BC = B. cenocepacia).The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.(EPS) S6 Fig. Model from no phage data, trained on only pairwise experimental data.Model fit predictions for 2-, 3-, and full 4-species community dynamics (solid lines) compared to experimental data (dashed lines) (PA = P. aeruginosa, SA = S. aureus, AB = A. baumannii, BC = B. Long time simulation of full community model shows shift in ecological outcomesgiven inclusion of phage.Simulation of the 4-species community gLV model over a long time scale reveals a qualitative shift in the outcome of the community when phage is present (PA = P. aeruginosa, SA = S. aureus, AB = A. baumannii, BC = B. cenocepacia).In the absence of phage (top), P. aeruginosa is the dominant competitor and only surviving species.In the presence of phage (bottom), the dominant competitor is eliminated, and we see competitive release of A. baumannii and S. aureus-maintaining 2 of the 3 non-targeted species in the community.Growth and interaction coefficients for simulation are from the model fits in Figs 9 and 10 and are shown for the wild-type PA14, no phage case (top left) in S8 Fig.For a detailed description of model parameterisation and simulation methods, see Methods and S1 Text.The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.(EPS) S8 Fig. Inferred interaction coefficients for the fitted gLV model describing full community dynamics, using 2-and 3-species experimental data with wild-type PA14 in the absence of phage (Fig 9) (PA = P. aeruginosa, SA = S. aureus, AB = A. baumannii, BC = B. cenocepacia).Heat map depicts β i,j coefficients (also labelled) scaled by P. aeruginosa intraspecific competition (β 1,1 = 1.2617×10 −8 , top left) corresponding to the wild-type PA14
The fold change difference of the individual community species not targeted by phage when comparing absolute densities in the presence of phage to the absence at the final experimental time point.Asterisks indicate higher final absolute density in the presence versus absence of phage (Wilcoxon signed rank exact test: * p < 0.05, ** p < 0.01, *** p < 0.001).The data and code required to generate this figure can be found in https://doi.org/10.6084/m9.figshare.24187284.